rm(list = ls())
library(ggplot2)
library(magrittr)
library(signal)
library(rootSolve)
devtools::load_all()

# 参数设置
beta <- .97
R <-  .02
WL <- 10

# 欧拉方程迭代收敛
cc <- B <- seq(-600,0,length.out = 80)
cc <- matrix(cc, ncol = 1)
for(i in 1:600){
  B_next <- (1+R) * B + WL - beta * (1+R) * cc[,ncol(cc)]
  c_next <- interp1(B,cc[,ncol(cc)], B_next,extrap = T)
  cc <- cbind(cc,c_next/(beta*(1+R)))
}
policy <- data.frame(B = B, cc = cc[,ncol(cc)])

# 按照政策函数迭代
cpath[1] <- interp1(policy$B, policy$cc, B[1])
for (i in 2:30) {
  cpath[i] <- beta * (1+R) * cpath[i-1]
}
data.frame(B = interp1(policy$cc,policy$B, cpath), cc = cpath)

# 按照边界条件迭代
